Reworking of Marcus’ code

Overview

  • A shortened version of the PT protocol is shown to have very bad internal consistency, test-retest reliability, convergent validity, and poor structural validity.
  • The issue is that it is a shortened version of the protocol. For this to be more widely accepted, we might need to run an additional study with a longer version of it.

TODO

  • remove PTTA from CFA? it was a separate exploratory hypothesis

Dependencies

library(readr)
library(dplyr)
library(tidyr)
library(stringr)
library(forcats)
library(irr)
library(psych)
library(lavaan)
library(janitor)
library(ggplot2)
library(ggrain)
library(scales)
library(mirt)
library(correlation)
library(see)
library(insight)
library(purrr)
library(knitr)
library(kableExtra)

options(scipen=999)

Load data

data_processed <- read_csv("../data/processed/data_processed.csv") 

data_processed_after_exclusions <- data_processed |>
  filter(exclude_master == "include")

Demographics

Sample sizes and exclusions

data_processed |>
  count(timepoint, exclude_master) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
timepoint exclude_master n
1 exclude 2
1 include 134
2 exclude 2
2 include 99

Analytic sample (timepoint 1)

data_processed_after_exclusions |>
  filter(timepoint == 1) |>
  count(gender) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
gender n
Non binary 1
female 91
male 42 
data_processed_after_exclusions |>
  filter(timepoint == 1) |>
  dplyr::summarize(mean_age = mean(age),
                   sd_age = sd(age),
                   min_age = min(age),
                   max_age = max(age)) |>
  pivot_longer(cols = everything()) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
name value
mean_age 38.23
sd_age 13.04
min_age 19.00
max_age 75.00 
data_processed_after_exclusions |>
  distinct(unique_id, .keep_all = TRUE) |>
  ggplot(aes(age)) +
  geom_histogram(binwidth = 5) +
  theme_linedraw()

Analyses

Test-retest

\(ICC_{2,1}\)

Results

data_reshaped <- data_processed_after_exclusions |>
  select(unique_id, timepoint, SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, PET_total, EYES_total, PTP_simple, PTP_reversed, PTP_doublereversed, PTP_total, PTTA_total) |>
  pivot_wider(id_cols = unique_id,
              names_from = timepoint,
              values_from = c(SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, PET_total, EYES_total, PTP_simple, PTP_reversed, PTP_doublereversed, PTP_total, PTTA_total))

tidy_icc_2_1 <- function(data, scale){
  res <- data |>
    select(starts_with(scale)) |>
    irr::icc(model = "twoway", type = "agreement", unit = "single",
             r0 = 0, conf.level = 0.95)
  
  tibble(scale = scale,
         icc_2.1 = res$value,
         ci_lower = res$lbound,
         ci_upper = res$ubound)
}

res <- bind_rows(
  tidy_icc_2_1(data_reshaped, "SITES"),
  tidy_icc_2_1(data_reshaped, "IRI_FS"),
  tidy_icc_2_1(data_reshaped, "IRI_EC"),
  tidy_icc_2_1(data_reshaped, "IRI_PT"),
  tidy_icc_2_1(data_reshaped, "IRI_PD"),
  tidy_icc_2_1(data_reshaped, "IRI_total"),
  tidy_icc_2_1(data_reshaped, "EQ_total"),
  tidy_icc_2_1(data_reshaped, "PET_total"),
  tidy_icc_2_1(data_reshaped, "EYES_total"),
  tidy_icc_2_1(data_reshaped, "PTP_simple"),
  tidy_icc_2_1(data_reshaped, "PTP_reversed"),
  tidy_icc_2_1(data_reshaped, "PTP_doublereversed"),
  tidy_icc_2_1(data_reshaped, "PTP_total"),
  tidy_icc_2_1(data_reshaped, "PTTA_total")
)
  
res |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
scale icc_2.1 ci_lower ci_upper
SITES 0.77 0.67 0.84
IRI_FS 0.88 0.82 0.92
IRI_EC 0.87 0.82 0.91
IRI_PT 0.81 0.73 0.87
IRI_PD 0.88 0.82 0.92
IRI_total 0.90 0.86 0.93
EQ_total 0.66 0.53 0.76
PET_total 0.86 0.80 0.91
EYES_total 0.73 0.63 0.81
PTP_simple 0.20 0.01 0.38
PTP_reversed 0.33 0.14 0.49
PTP_doublereversed 0.16 -0.04 0.35
PTP_total 0.08 -0.12 0.27
PTTA_total 0.69 0.54 0.80

Plots

round_half_up_trailing <- function(x, digits = 2) {
  x_rounded <- round_half_up(x, digits = digits)
  sprintf("%.*f", digits, x_rounded)
}

res_print <- res |>
  mutate(across(where(is.numeric), ~ round_half_up_trailing(.x, digits = 2))) |>
  mutate(string = paste0("= ", icc_2.1, ", 95% CI [", ci_lower, ", ", ci_upper, "]")) |>
  select(scale, string)

data_reshaped_plotting <- data_processed_after_exclusions |>
  select(unique_id, timepoint, SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, PET_total, EYES_total, PTP_total, PTTA_total) |>
  pivot_longer(cols = c(SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, PET_total, EYES_total, PTP_total, PTTA_total),
               names_to = "scale",
               values_to = "score")

data_reshaped_plotting |>
  filter(scale == "IRI_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Interpersonal Reactivity Index:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "IRI_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

data_reshaped_plotting |>
  filter(scale == "IRI_PT") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("IRI - Perspective Taking Subscale:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "IRI_PT") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

data_reshaped_plotting |>
  filter(scale == "SITES") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id", likert = TRUE) +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Single Item Trait Empathy Scale:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "SITES") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

data_reshaped_plotting |>
  filter(scale == "PET_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Pictoral Empathy Test:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "PET_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

data_reshaped_plotting |>
  filter(scale == "EQ_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Emotional Quotient:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "EQ_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

data_reshaped_plotting |>
  filter(scale == "EYES_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Reading the Mind in the Eyes Test:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "EYES_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

data_reshaped_plotting |>
  filter(scale == "PTTA_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Perspective Taking Task for Adults:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "PTTA_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

data_reshaped_plotting |>
  filter(scale == "PTP_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id", likert = TRUE) +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Perspective Taking Protocol:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "PTP_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

Internal consistency

Cronbach’s \(\alpha\)

data_subset <- data_processed_after_exclusions |>
  dplyr::select(timepoint, 
                starts_with("IRI_"), 
                starts_with("EQ_"),
                starts_with("PET_"),
                starts_with("EYES_"),
                starts_with("PTP_"),
                starts_with("PTTA_")) |>
  dplyr::select(-PTP_simple, -PTP_reversed, -PTP_doublereversed) |>
  dplyr::select(-IRI_FS, -IRI_EC, -IRI_PT, -IRI_PD) |>
  dplyr::select(!ends_with("_completeness") & !ends_with("_total"))

tidy_alpha <- function(data, time, scale){
  res <- data |>
    filter(timepoint == time) |>
    dplyr::select(starts_with(scale)) |>
    psych::alpha()
  
  res2 <- res$feldt 
  
  tibble(scale = scale,
         timepoint = time,
         alpha = res2$alpha$raw_alpha,
         ci_lower = res2$lower.ci$raw_alpha,
         ci_upper = res2$upper.ci$raw_alpha)
}

res <- bind_rows(
  tidy_alpha(data = data_subset, time = 1, scale = "IRI_"),
  tidy_alpha(data = data_subset, time = 2, scale = "IRI_"),
  tidy_alpha(data = data_subset, time = 1, scale = "EQ_"),
  tidy_alpha(data = data_subset, time = 2, scale = "EQ_"),
  tidy_alpha(data = data_subset, time = 1, scale = "PET_"),
  tidy_alpha(data = data_subset, time = 2, scale = "PET_"),
  tidy_alpha(data = data_subset, time = 1, scale = "EYES_"),
  tidy_alpha(data = data_subset, time = 2, scale = "EYES_"),
  tidy_alpha(data = data_subset, time = 1, scale = "PTP_"),
  tidy_alpha(data = data_subset, time = 2, scale = "PTP_"),
  tidy_alpha(data = data_subset, time = 1, scale = "PTTA_"),
  tidy_alpha(data = data_subset, time = 2, scale = "PTTA_")
)
## Some items ( IRI_19 IRI_24 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( EQ_4 EQ_8 EQ_10 EQ_12 EQ_14 EQ_15 EQ_18 EQ_21 EQ_27 EQ_28 EQ_29 EQ_32 EQ_34 EQ_46 EQ_48 EQ_49 EQ_50 EQ_57 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( EQ_4 EQ_8 EQ_10 EQ_12 EQ_14 EQ_15 EQ_18 EQ_21 EQ_27 EQ_28 EQ_29 EQ_32 EQ_46 EQ_48 EQ_49 EQ_50 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( EYES_31 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( PTP_9 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( PTP_1 PTP_8 PTP_9 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( PTTA_28 PTTA_29 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( PTTA_29 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
res |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
scale timepoint alpha ci_lower ci_upper
IRI_ 1 0.88 0.85 0.91
IRI_ 2 0.88 0.84 0.91
EQ_ 1 0.75 0.69 0.81
EQ_ 2 0.79 0.73 0.85
PET_ 1 0.92 0.89 0.94
PET_ 2 0.91 0.88 0.94
EYES_ 1 0.68 0.60 0.76
EYES_ 2 0.77 0.69 0.83
PTP_ 1 0.40 0.24 0.54
PTP_ 2 0.29 0.06 0.48
PTTA_ 1 0.87 0.83 0.90
PTTA_ 2 0.89 0.86 0.92
  • TODO: add subscales for IRI + other scales

Correlations

dat <- data_processed |>
  filter(timepoint == 1)

Scale correlations

scale_cors <- dat |>
  select(SITES, ends_with("_total")) |>
  rename_with(~ str_remove(.x, "_total")) |>
  correlation()

scale_cors |>
  as.data.frame() |>
  select(Parameter1, Parameter2, r, ci_low = CI_low, ci_high = CI_high, p) |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
Parameter1 Parameter2 r ci_low ci_high p
SITES IRI 0.57 0.44 0.67 p < .001
SITES EQ 0.31 0.15 0.45 p = 0.005
SITES PET 0.44 0.30 0.57 p < .001
SITES EYES 0.07 -0.10 0.24 p > .999
SITES PTP 0.03 -0.14 0.19 p > .999
SITES PTTA -0.05 -0.22 0.12 p > .999
IRI EQ 0.14 -0.03 0.30 p > .999
IRI PET 0.65 0.55 0.74 p < .001
IRI EYES 0.09 -0.08 0.25 p > .999
IRI PTP 0.05 -0.12 0.22 p > .999
IRI PTTA -0.16 -0.32 0.01 p = 0.837
EQ PET 0.11 -0.06 0.28 p > .999
EQ EYES -0.04 -0.21 0.13 p > .999
EQ PTP -0.06 -0.23 0.11 p > .999
EQ PTTA -0.11 -0.27 0.06 p > .999
PET EYES 0.11 -0.06 0.27 p > .999
PET PTP -0.01 -0.18 0.16 p > .999
PET PTTA -0.18 -0.34 -0.01 p = 0.579
EYES PTP 0.20 0.03 0.35 p = 0.352
EYES PTTA 0.18 0.01 0.34 p = 0.579
PTP PTTA 0.04 -0.13 0.21 p > .999 
scale_cors |>
  summary() |>
  plot(show_data = "points") +
  theme_linedraw() +
  ggtitle("")

Subscale correlations

subscale_cors <- dat |>
  select(SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, 
         PTP_simple, PTP_reversed, PTP_doublereversed, 
         ends_with("_total")) |>
  rename_with(~ str_remove(.x, "_total")) |>
  correlation() 

subscale_cors |>
  as.data.frame() |>
  select(Parameter1, Parameter2, r, ci_low = CI_low, ci_high = CI_high, p) |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
Parameter1 Parameter2 r ci_low ci_high p
SITES IRI_FS 0.33 0.17 0.48 p = 0.006
SITES IRI_EC 0.68 0.58 0.76 p < .001
SITES IRI_PT 0.54 0.41 0.65 p < .001
SITES IRI_PD 0.12 -0.05 0.28 p > .999
SITES PTP_simple 0.14 -0.03 0.30 p > .999
SITES PTP_reversed 0.02 -0.15 0.18 p > .999
SITES PTP_doublereversed -0.06 -0.23 0.11 p > .999
SITES IRI 0.57 0.44 0.67 p < .001
SITES EQ 0.31 0.15 0.45 p = 0.018
SITES PET 0.44 0.30 0.57 p < .001
SITES EYES 0.07 -0.10 0.24 p > .999
SITES PTP 0.03 -0.14 0.19 p > .999
SITES PTTA -0.05 -0.22 0.12 p > .999
IRI_FS IRI_EC 0.52 0.38 0.63 p < .001
IRI_FS IRI_PT 0.36 0.20 0.50 p = 0.002
IRI_FS IRI_PD 0.28 0.12 0.43 p = 0.059
IRI_FS PTP_simple 0.16 0.00 0.32 p > .999
IRI_FS PTP_reversed 0.11 -0.06 0.27 p > .999
IRI_FS PTP_doublereversed -0.03 -0.20 0.14 p > .999
IRI_FS IRI 0.79 0.71 0.84 p < .001
IRI_FS EQ 0.06 -0.11 0.22 p > .999
IRI_FS PET 0.50 0.36 0.62 p < .001
IRI_FS EYES 0.10 -0.07 0.27 p > .999
IRI_FS PTP 0.11 -0.06 0.28 p > .999
IRI_FS PTTA -0.09 -0.25 0.08 p > .999
IRI_EC IRI_PT 0.65 0.54 0.74 p < .001
IRI_EC IRI_PD 0.25 0.09 0.40 p = 0.211
IRI_EC PTP_simple 0.19 0.02 0.35 p > .999
IRI_EC PTP_reversed 0.04 -0.13 0.20 p > .999
IRI_EC PTP_doublereversed -0.04 -0.21 0.13 p > .999
IRI_EC IRI 0.84 0.78 0.88 p < .001
IRI_EC EQ 0.24 0.08 0.40 p = 0.287
IRI_EC PET 0.61 0.49 0.71 p < .001
IRI_EC EYES 0.05 -0.12 0.22 p > .999
IRI_EC PTP 0.07 -0.10 0.23 p > .999
IRI_EC PTTA -0.18 -0.34 -0.01 p > .999
IRI_PT IRI_PD -0.01 -0.18 0.16 p > .999
IRI_PT PTP_simple -0.01 -0.18 0.16 p > .999
IRI_PT PTP_reversed -0.06 -0.23 0.11 p > .999
IRI_PT PTP_doublereversed -0.06 -0.22 0.11 p > .999
IRI_PT IRI 0.66 0.55 0.75 p < .001
IRI_PT EQ 0.25 0.09 0.40 p = 0.211
IRI_PT PET 0.43 0.28 0.56 p < .001
IRI_PT EYES 0.08 -0.09 0.25 p > .999
IRI_PT PTP -0.08 -0.24 0.09 p > .999
IRI_PT PTTA -0.05 -0.22 0.12 p > .999
IRI_PD PTP_simple 0.04 -0.13 0.21 p > .999
IRI_PD PTP_reversed 0.01 -0.16 0.18 p > .999
IRI_PD PTP_doublereversed -0.01 -0.18 0.16 p > .999
IRI_PD IRI 0.56 0.43 0.66 p < .001
IRI_PD EQ -0.13 -0.29 0.04 p > .999
IRI_PD PET 0.32 0.16 0.47 p = 0.010
IRI_PD EYES 0.01 -0.16 0.18 p > .999
IRI_PD PTP 0.02 -0.15 0.19 p > .999
IRI_PD PTTA -0.14 -0.30 0.03 p > .999
PTP_simple PTP_reversed 0.32 0.16 0.46 p = 0.012
PTP_simple PTP_doublereversed 0.05 -0.12 0.22 p > .999
PTP_simple IRI 0.15 -0.02 0.31 p > .999
PTP_simple EQ 0.06 -0.11 0.23 p > .999
PTP_simple PET 0.23 0.06 0.38 p = 0.465
PTP_simple EYES 0.26 0.10 0.41 p = 0.145
PTP_simple PTP 0.59 0.47 0.69 p < .001
PTP_simple PTTA 0.06 -0.11 0.22 p > .999
PTP_reversed PTP_doublereversed -0.06 -0.23 0.11 p > .999
PTP_reversed IRI 0.04 -0.13 0.21 p > .999
PTP_reversed EQ -0.11 -0.27 0.06 p > .999
PTP_reversed PET -0.10 -0.27 0.07 p > .999
PTP_reversed EYES 0.19 0.02 0.35 p > .999
PTP_reversed PTP 0.77 0.70 0.83 p < .001
PTP_reversed PTTA 0.06 -0.11 0.23 p > .999
PTP_doublereversed IRI -0.05 -0.21 0.12 p > .999
PTP_doublereversed EQ 0.00 -0.17 0.17 p > .999
PTP_doublereversed PET -0.02 -0.19 0.15 p > .999
PTP_doublereversed EYES -0.04 -0.21 0.13 p > .999
PTP_doublereversed PTP 0.50 0.36 0.62 p < .001
PTP_doublereversed PTTA -0.04 -0.21 0.13 p > .999
IRI EQ 0.14 -0.03 0.30 p > .999
IRI PET 0.65 0.55 0.74 p < .001
IRI EYES 0.09 -0.08 0.25 p > .999
IRI PTP 0.05 -0.12 0.22 p > .999
IRI PTTA -0.16 -0.32 0.01 p > .999
EQ PET 0.11 -0.06 0.28 p > .999
EQ EYES -0.04 -0.21 0.13 p > .999
EQ PTP -0.06 -0.23 0.11 p > .999
EQ PTTA -0.11 -0.27 0.06 p > .999
PET EYES 0.11 -0.06 0.27 p > .999
PET PTP -0.01 -0.18 0.16 p > .999
PET PTTA -0.18 -0.34 -0.01 p > .999
EYES PTP 0.20 0.03 0.35 p > .999
EYES PTTA 0.18 0.01 0.34 p > .999
PTP PTTA 0.04 -0.13 0.21 p > .999 
# code repeated in order to use the rename, which screws up the table but is needed for the plot
dat |>
  select(SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, 
         PTP_simple, PTP_reversed, PTP_doublereversed, 
         ends_with("_total")) |>
  rename_with(~ str_remove(.x, "_total")) |>
  rename_with(~ str_replace(.x, "_", "\n")) |>
  correlation()  |>
  summary() |>
  plot(show_data = "points") +
  theme_linedraw() +
  ggtitle("")

SEM

dat <- data_processed |>
  filter(timepoint == 1) |>
  select(SITES, ends_with("_total"), IRI_PT)

# model <- '
#          PT =~ SITES + IRI_total + EQ_total + PET_total + EYES_total + PTP_total + PTTA_total
#          '

model <- '
         PT =~ SITES + IRI_total + EQ_total + PET_total + EYES_total + PTP_total
         '

fit <- cfa(model = model, 
           data = dat)

fitMeasures(fit, fit.measures = c("cfi", "tli", "rmsea", "srmr", "chisq", "df", "pvalue"))
##    cfi    tli  rmsea   srmr  chisq     df pvalue 
##  0.976  0.961  0.051  0.048 12.145  9.000  0.205
standardizedSolution(fit) |>
  filter(op != "~~") |>
  select(scale = rhs, est.std, ci.lower, ci.upper, p = pvalue) |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
scale est.std ci.lower ci.upper p
SITES 0.64 0.52 0.76 p < .001
IRI_total 0.89 0.79 1.00 p < .001
EQ_total 0.19 0.02 0.37 p = 0.033
PET_total 0.73 0.61 0.84 p < .001
EYES_total 0.11 -0.07 0.29 p = 0.245
PTP_total 0.04 -0.14 0.22 p = 0.678
  • The overall model fits adequately, suggesting that the scales are collectively measuring a latent factor (empathy/PT).
  • Non significant loading of PTP onto the latent factor suggests it is not a good measure of this latent factor of empathy/PT.

Rasch analysis of PTP

Unidimensional

# subset data
dat <- data_processed_after_exclusions |>
  filter(timepoint == 1) |>
  select(starts_with("PTP_")) |>
  select(-PTP_simple, -PTP_reversed, -PTP_doublereversed, -PTP_completeness, -PTP_total)

# fit model
#uni_model <- mirt(dat, model = 1, itemtype = "2PL")
uni_model <- mirt(dat, model = 1, itemtype = "Rasch")
## Iteration: 1, Log-Lik: -507.838, Max-Change: 0.26044Iteration: 2, Log-Lik: -507.188, Max-Change: 0.05781Iteration: 3, Log-Lik: -506.923, Max-Change: 0.04393Iteration: 4, Log-Lik: -506.754, Max-Change: 0.03920Iteration: 5, Log-Lik: -506.607, Max-Change: 0.02715Iteration: 6, Log-Lik: -506.526, Max-Change: 0.02191Iteration: 7, Log-Lik: -506.477, Max-Change: 0.02233Iteration: 8, Log-Lik: -506.420, Max-Change: 0.01476Iteration: 9, Log-Lik: -506.391, Max-Change: 0.01222Iteration: 10, Log-Lik: -506.372, Max-Change: 0.01159Iteration: 11, Log-Lik: -506.353, Max-Change: 0.00862Iteration: 12, Log-Lik: -506.342, Max-Change: 0.00730Iteration: 13, Log-Lik: -506.336, Max-Change: 0.01047Iteration: 14, Log-Lik: -506.326, Max-Change: 0.00530Iteration: 15, Log-Lik: -506.322, Max-Change: 0.00450Iteration: 16, Log-Lik: -506.319, Max-Change: 0.00427Iteration: 17, Log-Lik: -506.316, Max-Change: 0.00333Iteration: 18, Log-Lik: -506.314, Max-Change: 0.00288Iteration: 19, Log-Lik: -506.313, Max-Change: 0.00417Iteration: 20, Log-Lik: -506.311, Max-Change: 0.00216Iteration: 21, Log-Lik: -506.311, Max-Change: 0.00185Iteration: 22, Log-Lik: -506.310, Max-Change: 0.00174Iteration: 23, Log-Lik: -506.310, Max-Change: 0.00139Iteration: 24, Log-Lik: -506.309, Max-Change: 0.00121Iteration: 25, Log-Lik: -506.309, Max-Change: 0.00195Iteration: 26, Log-Lik: -506.309, Max-Change: 0.00092Iteration: 27, Log-Lik: -506.309, Max-Change: 0.00080Iteration: 28, Log-Lik: -506.309, Max-Change: 0.00075Iteration: 29, Log-Lik: -506.308, Max-Change: 0.00060Iteration: 30, Log-Lik: -506.308, Max-Change: 0.00053Iteration: 31, Log-Lik: -506.308, Max-Change: 0.00093Iteration: 32, Log-Lik: -506.308, Max-Change: 0.00040Iteration: 33, Log-Lik: -506.308, Max-Change: 0.00035Iteration: 34, Log-Lik: -506.308, Max-Change: 0.00032Iteration: 35, Log-Lik: -506.308, Max-Change: 0.00027Iteration: 36, Log-Lik: -506.308, Max-Change: 0.00024Iteration: 37, Log-Lik: -506.308, Max-Change: 0.00032Iteration: 38, Log-Lik: -506.308, Max-Change: 0.00018Iteration: 39, Log-Lik: -506.308, Max-Change: 0.00016Iteration: 40, Log-Lik: -506.308, Max-Change: 0.00017Iteration: 41, Log-Lik: -506.308, Max-Change: 0.00012Iteration: 42, Log-Lik: -506.308, Max-Change: 0.00011Iteration: 43, Log-Lik: -506.308, Max-Change: 0.00009
# model fit stats
M2(uni_model) |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
M2 df p RMSEA RMSEA_5 RMSEA_95 SRMSR TLI CFI
stats 99.14 35 p < .001 0.12 0.09 0.14 0.15 0.37 0.39 
# item parameters
coefs_uni <- coef(uni_model, IRTpars = TRUE) # Extract item parameters in IRT metric

item_parameters <- coefs_uni  |>
  discard(names(coefs_uni) == "GroupPars") %>% # Exclude GroupPars if needed
  imap_dfr(~ as_tibble(.x, rownames = "parameter") %>% 
             mutate(PTP = .y)) %>%
  relocate(PTP) |>
  rename(item = PTP) |>
  mutate(subscale = case_when(item %in% c("PTP_1", "PTP_2", "PTP_3") ~ "simple",
                              item %in% c("PTP_4", "PTP_5", "PTP_6") ~ "reversed",
                              item %in% c("PTP_7", "PTP_8", "PTP_9") ~ "doublereversed"),
         subscale = fct_relevel(subscale, "simple", "reversed", "doublereversed")) |>
  select(item, subscale, difficulty = b)

item_parameters |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
item subscale difficulty
PTP_1 simple -2.56
PTP_2 simple -3.53
PTP_3 simple -4.06
PTP_4 reversed -2.07
PTP_5 reversed -1.70
PTP_6 reversed -1.76
PTP_7 doublereversed -1.04
PTP_8 doublereversed 0.38
PTP_9 doublereversed 0.38 
item_parameters |>
  group_by(subscale) |>
  dplyr::summarize(mean_difficulty = mean(difficulty)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
subscale mean_difficulty
simple -3.38
reversed -1.84
doublereversed -0.10
  • p < 0.05 indicates potential model misfit.
  • RMSEA < 0.05: Good fit; RMSEA between 0.05 and 0.08: Acceptable fit; RMSEA > 0.10: Poor fit.

Multidimensional

# Model specification: 'F1' loads on items 1-3, 'F2' on items 4-6, 'F3' on items 7-9.
# three_factor_model <- '
#                       F1 = 1-3
#                       F2 = 4-6
#                       F3 = 7-9
#                       '

# PTP_7 was incorrectly implemented as a single reversed trial in the survey instead of double. Corrected here.
three_factor_model <- '
                      F1 = 1-3
                      F2 = 4-7
                      F3 = 8-9
                      '

# Fit the multidimensional model
multi_model <- mirt(dat, model = three_factor_model, itemtype = "Rasch")
## Iteration: 1, Log-Lik: -495.722, Max-Change: 0.29425Iteration: 2, Log-Lik: -495.111, Max-Change: 0.06521Iteration: 3, Log-Lik: -494.759, Max-Change: 0.06460Iteration: 4, Log-Lik: -494.437, Max-Change: 0.07189Iteration: 5, Log-Lik: -494.103, Max-Change: 0.06468Iteration: 6, Log-Lik: -493.839, Max-Change: 0.06079Iteration: 7, Log-Lik: -493.631, Max-Change: 0.07020Iteration: 8, Log-Lik: -493.339, Max-Change: 0.05848Iteration: 9, Log-Lik: -493.143, Max-Change: 0.05254Iteration: 10, Log-Lik: -492.970, Max-Change: 0.05473Iteration: 11, Log-Lik: -492.786, Max-Change: 0.04652Iteration: 12, Log-Lik: -492.641, Max-Change: 0.04478Iteration: 13, Log-Lik: -492.531, Max-Change: 0.05879Iteration: 14, Log-Lik: -492.350, Max-Change: 0.04940Iteration: 15, Log-Lik: -492.234, Max-Change: 0.04892Iteration: 16, Log-Lik: -492.127, Max-Change: 0.05933Iteration: 17, Log-Lik: -492.008, Max-Change: 0.05264Iteration: 18, Log-Lik: -491.913, Max-Change: 0.05235Iteration: 19, Log-Lik: -491.845, Max-Change: 0.07219Iteration: 20, Log-Lik: -491.708, Max-Change: 0.05776Iteration: 21, Log-Lik: -491.625, Max-Change: 0.05584Iteration: 22, Log-Lik: -491.547, Max-Change: 0.06859Iteration: 23, Log-Lik: -491.458, Max-Change: 0.05926Iteration: 24, Log-Lik: -491.387, Max-Change: 0.05795Iteration: 25, Log-Lik: -491.339, Max-Change: 0.08389Iteration: 26, Log-Lik: -491.227, Max-Change: 0.06386Iteration: 27, Log-Lik: -491.164, Max-Change: 0.06009Iteration: 28, Log-Lik: -491.104, Max-Change: 0.07454Iteration: 29, Log-Lik: -491.035, Max-Change: 0.06260Iteration: 30, Log-Lik: -490.982, Max-Change: 0.06005Iteration: 31, Log-Lik: -490.948, Max-Change: 0.09058Iteration: 32, Log-Lik: -490.858, Max-Change: 0.06580Iteration: 33, Log-Lik: -490.812, Max-Change: 0.06003Iteration: 34, Log-Lik: -490.767, Max-Change: 0.07423Iteration: 35, Log-Lik: -490.717, Max-Change: 0.06085Iteration: 36, Log-Lik: -490.679, Max-Change: 0.05712Iteration: 37, Log-Lik: -490.655, Max-Change: 0.08918Iteration: 38, Log-Lik: -490.588, Max-Change: 0.06180Iteration: 39, Log-Lik: -490.555, Max-Change: 0.05433Iteration: 40, Log-Lik: -490.524, Max-Change: 0.06624Iteration: 41, Log-Lik: -490.490, Max-Change: 0.05310Iteration: 42, Log-Lik: -490.464, Max-Change: 0.04881Iteration: 43, Log-Lik: -490.447, Max-Change: 0.07810Iteration: 44, Log-Lik: -490.402, Max-Change: 0.05169Iteration: 45, Log-Lik: -490.380, Max-Change: 0.04397Iteration: 46, Log-Lik: -490.359, Max-Change: 0.05221Iteration: 47, Log-Lik: -490.338, Max-Change: 0.04122Iteration: 48, Log-Lik: -490.321, Max-Change: 0.03719Iteration: 49, Log-Lik: -490.310, Max-Change: 0.06063Iteration: 50, Log-Lik: -490.281, Max-Change: 0.03854Iteration: 51, Log-Lik: -490.267, Max-Change: 0.03170Iteration: 52, Log-Lik: -490.255, Max-Change: 0.03618Iteration: 53, Log-Lik: -490.242, Max-Change: 0.02843Iteration: 54, Log-Lik: -490.232, Max-Change: 0.02536Iteration: 55, Log-Lik: -490.224, Max-Change: 0.04191Iteration: 56, Log-Lik: -490.208, Max-Change: 0.02590Iteration: 57, Log-Lik: -490.199, Max-Change: 0.02072Iteration: 58, Log-Lik: -490.192, Max-Change: 0.02300Iteration: 59, Log-Lik: -490.184, Max-Change: 0.01804Iteration: 60, Log-Lik: -490.178, Max-Change: 0.01596Iteration: 61, Log-Lik: -490.174, Max-Change: 0.02666Iteration: 62, Log-Lik: -490.164, Max-Change: 0.01618Iteration: 63, Log-Lik: -490.159, Max-Change: 0.01261Iteration: 64, Log-Lik: -490.155, Max-Change: 0.01377Iteration: 65, Log-Lik: -490.151, Max-Change: 0.01071Iteration: 66, Log-Lik: -490.147, Max-Change: 0.00943Iteration: 67, Log-Lik: -490.145, Max-Change: 0.01601Iteration: 68, Log-Lik: -490.140, Max-Change: 0.00952Iteration: 69, Log-Lik: -490.137, Max-Change: 0.00736Iteration: 70, Log-Lik: -490.134, Max-Change: 0.00782Iteration: 71, Log-Lik: -490.132, Max-Change: 0.00619Iteration: 72, Log-Lik: -490.130, Max-Change: 0.00542Iteration: 73, Log-Lik: -490.128, Max-Change: 0.00920Iteration: 74, Log-Lik: -490.126, Max-Change: 0.00544Iteration: 75, Log-Lik: -490.124, Max-Change: 0.00415Iteration: 76, Log-Lik: -490.123, Max-Change: 0.00441Iteration: 77, Log-Lik: -490.121, Max-Change: 0.00347Iteration: 78, Log-Lik: -490.120, Max-Change: 0.00304Iteration: 79, Log-Lik: -490.120, Max-Change: 0.00528Iteration: 80, Log-Lik: -490.118, Max-Change: 0.00307Iteration: 81, Log-Lik: -490.117, Max-Change: 0.00234Iteration: 82, Log-Lik: -490.116, Max-Change: 0.00253Iteration: 83, Log-Lik: -490.116, Max-Change: 0.00199Iteration: 84, Log-Lik: -490.115, Max-Change: 0.00169Iteration: 85, Log-Lik: -490.115, Max-Change: 0.00283Iteration: 86, Log-Lik: -490.114, Max-Change: 0.00185Iteration: 87, Log-Lik: -490.113, Max-Change: 0.00126Iteration: 88, Log-Lik: -490.113, Max-Change: 0.00127Iteration: 89, Log-Lik: -490.112, Max-Change: 0.00104Iteration: 90, Log-Lik: -490.112, Max-Change: 0.00096Iteration: 91, Log-Lik: -490.112, Max-Change: 0.00157Iteration: 92, Log-Lik: -490.111, Max-Change: 0.00092Iteration: 93, Log-Lik: -490.111, Max-Change: 0.00072Iteration: 94, Log-Lik: -490.111, Max-Change: 0.00088Iteration: 95, Log-Lik: -490.111, Max-Change: 0.00056Iteration: 96, Log-Lik: -490.110, Max-Change: 0.00052Iteration: 97, Log-Lik: -490.110, Max-Change: 0.00074Iteration: 98, Log-Lik: -490.110, Max-Change: 0.00049Iteration: 99, Log-Lik: -490.110, Max-Change: 0.00046Iteration: 100, Log-Lik: -490.110, Max-Change: 0.00046Iteration: 101, Log-Lik: -490.110, Max-Change: 0.00040Iteration: 102, Log-Lik: -490.110, Max-Change: 0.00038Iteration: 103, Log-Lik: -490.110, Max-Change: 0.00037Iteration: 104, Log-Lik: -490.109, Max-Change: 0.00033Iteration: 105, Log-Lik: -490.109, Max-Change: 0.00031Iteration: 106, Log-Lik: -490.109, Max-Change: 0.00029Iteration: 107, Log-Lik: -490.109, Max-Change: 0.00029Iteration: 108, Log-Lik: -490.109, Max-Change: 0.00027Iteration: 109, Log-Lik: -490.109, Max-Change: 0.00027Iteration: 110, Log-Lik: -490.109, Max-Change: 0.00025Iteration: 111, Log-Lik: -490.109, Max-Change: 0.00023Iteration: 112, Log-Lik: -490.109, Max-Change: 0.00021Iteration: 113, Log-Lik: -490.109, Max-Change: 0.00021Iteration: 114, Log-Lik: -490.109, Max-Change: 0.00019Iteration: 115, Log-Lik: -490.109, Max-Change: 0.00019Iteration: 116, Log-Lik: -490.109, Max-Change: 0.00018Iteration: 117, Log-Lik: -490.109, Max-Change: 0.00017Iteration: 118, Log-Lik: -490.109, Max-Change: 0.00015Iteration: 119, Log-Lik: -490.109, Max-Change: 0.00015Iteration: 120, Log-Lik: -490.109, Max-Change: 0.00014Iteration: 121, Log-Lik: -490.109, Max-Change: 0.00014Iteration: 122, Log-Lik: -490.109, Max-Change: 0.00013Iteration: 123, Log-Lik: -490.109, Max-Change: 0.00012Iteration: 124, Log-Lik: -490.109, Max-Change: 0.00012Iteration: 125, Log-Lik: -490.109, Max-Change: 0.00012Iteration: 126, Log-Lik: -490.109, Max-Change: 0.00011Iteration: 127, Log-Lik: -490.109, Max-Change: 0.00010Iteration: 128, Log-Lik: -490.109, Max-Change: 0.00010
# Compare models using anova 
anova(uni_model, multi_model)
##                  AIC    SABIC       HQ      BIC   logLik     X2 df p
## uni_model   1032.616 1029.962 1044.392 1061.595 -506.308            
## multi_model 1004.217 1001.032 1018.348 1038.991 -490.109 32.399  2 0
# > multi provides better fit

# model fit stats
M2(multi_model) |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 3))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
M2 df p RMSEA RMSEA_5 RMSEA_95 SRMSR TLI CFI
stats 55.659 33 p = 0.008 0.072 0.037 0.103 0.12 0.764 0.784 
# item parameters
coefs_multi <- coef(multi_model, IRTpars = TRUE) # Extract item parameters in IRT metric

item_parameters_multi <- coefs_multi  |>
  discard(names(coefs_multi) == "GroupPars") %>% # Exclude GroupPars if needed
  imap_dfr(~ as_tibble(.x, rownames = "parameter") %>% 
             mutate(PTP = .y)) %>%
  relocate(PTP) |>
  rename(item = PTP) |>
  mutate(subscale = case_when(item %in% c("PTP_1", "PTP_2", "PTP_3") ~ "simple",
                              item %in% c("PTP_4", "PTP_5", "PTP_6") ~ "reversed",
                              item %in% c("PTP_7", "PTP_8", "PTP_9") ~ "doublereversed"),
         subscale = fct_relevel(subscale, "simple", "reversed", "doublereversed")) |>
  select(item, subscale, difficulty = b)

item_parameters_multi |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
item subscale difficulty
PTP_1 simple -3.66
PTP_2 simple -4.91
PTP_3 simple -5.55
PTP_4 reversed -2.50
PTP_5 reversed -2.07
PTP_6 reversed -2.13
PTP_7 doublereversed -1.28
PTP_8 doublereversed 0.45
PTP_9 doublereversed 0.45 
item_parameters_multi |>
  group_by(subscale) |>
  dplyr::summarize(mean_difficulty = mean(difficulty)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
subscale mean_difficulty
simple -4.71
reversed -2.23
doublereversed -0.12
  • p < 0.05 indicates potential model misfit.
  • RMSEA < 0.05: Good fit; RMSEA between 0.05 and 0.08: Acceptable fit; RMSEA > 0.10: Poor fit.

Session info

sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 15.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Zurich
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] kableExtra_1.4.0  knitr_1.49        purrr_1.0.2       insight_0.19.10  
##  [5] see_0.8.3         correlation_0.8.4 mirt_1.40         lattice_0.22-6   
##  [9] scales_1.3.0      ggrain_0.0.4      ggplot2_3.5.1     janitor_2.2.0    
## [13] lavaan_0.6-17     psych_2.4.6.26    irr_0.84.1        lpSolve_5.6.20   
## [17] forcats_1.0.0     stringr_1.5.1     tidyr_1.3.1       dplyr_1.1.4      
## [21] readr_2.1.5      
## 
## loaded via a namespace (and not attached):
##  [1] mnormt_2.1.1         pbapply_1.7-2        polynom_1.4-1       
##  [4] gridExtra_2.3        permute_0.9-7        sandwich_3.1-0      
##  [7] rlang_1.1.4          magrittr_2.0.3       multcomp_1.4-25     
## [10] snakecase_0.11.1     compiler_4.3.3       mgcv_1.9-1          
## [13] systemfonts_1.0.6    vctrs_0.6.5          quadprog_1.5-8      
## [16] pkgconfig_2.0.3      crayon_1.5.2         fastmap_1.2.0       
## [19] labeling_0.4.3       pbivnorm_0.6.0       utf8_1.2.4          
## [22] rmarkdown_2.29       tzdb_0.4.0           bit_4.0.5           
## [25] xfun_0.49            cachem_1.1.0         jsonlite_1.8.9      
## [28] gghalves_0.1.4       Deriv_4.1.3          parallel_4.3.3      
## [31] cluster_2.1.6        R6_2.5.1             bslib_0.8.0         
## [34] stringi_1.8.4        lubridate_1.9.4      jquerylib_0.1.4     
## [37] estimability_1.5.1   Rcpp_1.0.12          zoo_1.8-12          
## [40] parameters_0.21.6    Matrix_1.6-5         splines_4.3.3       
## [43] timechange_0.3.0     tidyselect_1.2.1     rstudioapi_0.17.1   
## [46] yaml_2.3.10          vegan_2.6-4          codetools_0.2-20    
## [49] dcurver_0.9.2        tibble_3.2.1         withr_3.0.2         
## [52] bayestestR_0.13.2    coda_0.19-4.1        evaluate_1.0.1      
## [55] survival_3.5-8       ggpp_0.5.6           xml2_1.3.6          
## [58] pillar_1.9.0         generics_0.1.3       vroom_1.6.5         
## [61] hms_1.1.3            munsell_0.5.1        xtable_1.8-4        
## [64] glue_1.8.0           emmeans_1.10.2       tools_4.3.3         
## [67] mvtnorm_1.2-4        grid_4.3.3           datawizard_0.10.0   
## [70] colorspace_2.1-1     nlme_3.1-164         cli_3.6.3           
## [73] fansi_1.0.6          viridisLite_0.4.2    svglite_2.1.3       
## [76] gtable_0.3.6         sass_0.4.9           digest_0.6.37       
## [79] GPArotation_2023.8-1 TH.data_1.1-2        farver_2.1.2        
## [82] htmltools_0.5.8.1    lifecycle_1.0.4      bit64_4.0.5         
## [85] MASS_7.3-60.0.1